# include "scene.h"
# include "discr.h"
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <errno.h>
# include <math.h>


# ifdef SOLVE_ENERGY
extern double temp_advance(struct grid *g, const struct arr_struct *depVar);
# endif

extern void init_cond(struct grid *g);
extern short print_result(unsigned int t_cnt, const struct grid *g);
extern unsigned int vel_advance(struct grid *g, unsigned int *pcnt_max, struct arr_struct *depVar, struct arr_struct (*memFace)[3]);
extern double inlet_velocity(double r, double t, short index);
extern double dens(double T);
extern double kvisc(double T);
extern double thm_cond(double T);
extern short print_recov(const struct grid *g);
extern double orlansky_outlet(struct grid *g, unsigned int t_cnt, const struct arr_struct *memOut);
extern FILE *file_open(FILE *lptr, const char *file, unsigned int t_cnt, unsigned int nfreq);

static double heat_balance(const struct grid *g, double E_absr);


extern short solve(struct grid *g, struct arr_struct *depVar, struct arr_struct (*memFace)[3], struct arr_struct *memOut)
 {
  unsigned int t_cnt, pcnt_max, pcnt_total;
  short sst_cnt, div_cnt;
  double pm_accum, pE_accum, E_absr;
  FILE *tptr;

# ifdef SOLVE_ENERGY
  g->vhsrc_flag = VHSRC_FLAG;
# endif

  g->density0 = dens(T_INLET);
  g->kvisc0 = kvisc(T_INLET);
  g->Re0 = 0.5*W_MAX*d_NOZZLE/g->kvisc0;
  fprintf(g->lptr, "\nReynolds Number at Inlet, Re0 = %E\n", g->Re0);
  g->thmcond0 = thm_cond(T_INLET);

  init_cond(g);

/* =================== Calculation at different time-steps ==================== */

  pcnt_total = sst_cnt = div_cnt = 0;
  tptr = (FILE *)0;
  t_cnt = 0;

  while(g->t <= tMAX)
    {
#   ifdef ittr_RES_FILE
     if(t_cnt%IRFILE_FREQ == 0)
          {
           if((tptr != (FILE *)0) && (tptr != stdout)) fclose(tptr);
           tptr = file_open(g->lptr, ittr_RES_FILE, t_cnt, IRFILE_FREQ);
           setlinebuf(tptr);
          }
#    endif

  /* Print in the itteration file per ittr_ENTRY_FREQ entries.  */
     g->tptr = (((t_cnt < ittr_PRINTALL_CNT) || ((t_cnt+1)%ittr_ENTRY_FREQ == 0)) ? tptr:(FILE *)0);

     if(g->tptr != (FILE *)0) fprintf(g->tptr, "\n\n<=======[[ t_cnt = %u, t = %E Sec ]]=======>\n", t_cnt, g->t);

     pm_accum = orlansky_outlet(g, t_cnt, memOut);
     if(g->tptr != (FILE *)0) fprintf(g->tptr, "\tpm_accum = %E%c\n", pm_accum, '%');

     pcnt_total += vel_advance(g, &pcnt_max, depVar, memFace);

#    ifdef SOLVE_ENERGY
     E_absr = temp_advance(g, depVar);
     pE_accum = heat_balance((const struct grid *)g, E_absr);
     if(g->tptr != (FILE *)0) fprintf(g->tptr, "\tpE_accum = %E%c\n", pE_accum, '%');
#    endif

     if(pcnt_max == POISS_CNT_MIN) ++sst_cnt; else sst_cnt = 0;
     if(pcnt_max == POISS_CNT_MAX+1) ++div_cnt; else div_cnt = 0;

     t_cnt++;
     g->t += g->dt;

  /* Printing result and recovery files. */
     if( t_cnt%RESF_FREQ == 0 || sst_cnt == SST_CNT_MAX ) print_result(t_cnt, g);

     if(pcnt_total >= RECOVF_POISS_CNT) if(pcnt_max <= POISS_CNT_MAX)
          {
           print_recov(g);
           pcnt_total = 0;
          }

     if(sst_cnt == SST_CNT_MAX)
           {
            if(g->tptr != (FILE *)0) fprintf(g->tptr, "(-: STEADY-STATE solution! :-)\n");
            fprintf(g->lptr, "(-: STEADY-STATE solution! :-)\n");
            break;
           }

     if(div_cnt == DIV_CNT_MAX)
           {
            if(g->tptr != (FILE *)0) fprintf(g->tptr, "***DIVERGED solution***\n");
            fprintf(g->lptr, "***DIVERGED Solution***\n");
            break;
           }
    }
/* --------------------------------------Time-loop ends------------------------------------------ */

  if((tptr != (FILE *)0) && (tptr != stdout)) fclose(tptr);

  return g->t > tMAX ? 1:-1;
 }


# ifdef SOLVE_ENERGY
static double heat_balance(const struct grid *g, double E_absr)
    {
     double E_inlet, E_outlet, E_vhsrc, E_heater, *th, *r/*, ***T*/;
     unsigned short i, j, ii, jj, nth/*, nr*/;

     nth = g->nth; th = g->th; /*nr = g->nr[2]; */r = g->r[2];/* T = g->T[2];*/

  /* Determining total energy created by the volume-heat-source. */
     E_vhsrc = ( g->vhsrc_flag == 'y' ? (th[NthG+nth]-th[NthG])*(\
          d0_CYL*d0_CYL*h0_CYL + 4*g->r[2][NrG+g->nr[2]]*g->r[2][NrG+g->nr[2]]*h1_CYL\
          + 4*(g->r[2][NrG+g->nr[2]]*g->r[2][NrG+g->nr[2]]-g->r[2][NrG]*g->r[2][NrG])*h2_CYL\
                             )*g->density0*Cp_WATER*(EXPECTED_T_RISE/EXPECTED_TIME)*g->dt/8:0 );

  /* Calculating heat input because of heater on the top surface. */
     E_heater = g->dt*TOP_HEATER_HF*0.5*g->r_heater*g->r_heater*(th[NthG+g->nth]-th[NthG]);

  /* Calculating the energy coming into the domain through the inlet in time t~t+dt. */
     E_inlet = 0.25*Cp_WATER*T_INLET*g->density0*r_NOZZLE*r_NOZZLE*W_MAX\
       *(g->dt-tFRAC*(exp(log(1-wFRAC)*(g->t+g->dt)/tFRAC) - exp(log(1-wFRAC)*g->t/tFRAC))/log(1-wFRAC))*(th[NthG+g->nth]-th[NthG]);

  /* Calculating the heat gone out of the domain through outlet between t~t+dt. */
/*     for(E_outlet = 0, i = 0, ii = NrG+nr-g->nr_gap; i < g->nr_gap; i++, ii++) for(j = 0, jj = NthG; j < nth; j++, jj++)
     E_outlet += -0.0625*(r[ii+1]+r[ii])*(r[ii+1]-r[ii])*(th[jj+1]-th[jj])*(g->outBC.w[0][i][j][0]+g->outBC.tnext_wg[i][j])\
                     *(g->outBC.T[0][i][j][0]+g->outBC.T[0][i][j][1]+g->outBC.tnext_Tg[i][j]+g->T[2][ii][jj][NzG]); */

     for(E_outlet = 0, i = 0; i < g->nr[2]; i++) for(ii = NrG+i, j = 0, jj = NthG; j < nth; j++, jj++)
     E_outlet += -0.0625*(r[ii+1]+r[ii])*(r[ii+1]-r[ii])*(th[jj+1]-th[jj])*(g->outBC.w[0][i][j][0]+g->outBC.tnext_wg[i][j])\
                     *(g->outBC.T[0][i][j][0]+g->outBC.T[0][i][j][1]+g->outBC.tnext_Tg[i][j]+g->T[2][ii][jj][NzG]);

     E_outlet *= g->density0*Cp_WATER*g->dt;

  /* Return % heat-accumulation in the domain between time t~t+dt. */
     return (1-(E_outlet+E_absr)/(E_inlet+E_heater+E_vhsrc))*100;
    }
# endif
